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I. INTRODUCTION. 



This report presents two great circle navigation algorithms, a closest point of approach 
(CPA) algorithm and an intercept algorithm. It also presents an implementation program 
that is written in the BASIC language for an IBM PC. Instead of using classical spherical 
geometry or the general spherical triangle [Ref. 1], these algorithms incorporate rectangular 
coordinates and vectors on the surface of the sphere. Portions of the vector formalism were 
suggested by Reference 2. 

The intent of the report is to provide algorithms for spherical earth models that can 
be used for situations in which flat earth models are not appropriate, but that do not 
require the sophistication of rotating earth models. 



H. RECTANGULAR COORDINATES AND VECTORS ON A SPHERE 

In a spherical earth model, a point P on the earth’s surface can be located by a position 
vector p = xi+yj+zk in a rectangular coordinate system with origin at the earth’s center. 
In matrix form, and with x, y and z expressed in spherical coordinates, 



P = 



r cos <f> cos A 
rcos^sinA 

ra\n<f> 



where <f> is the latitude and A is the longitude at the point and r is the distance of the 
point from the earth’s center. See Figure 1. 

In terms of the unit vector x where 



(1) 



p = rx. We can think of x as the unit vector normal to the surface at the point. It is 
convenient to define two other unit vectors in the tangent plane to the earth’s surface at 
the point. These are local north n and local east e, defined by 





'xj ' 




’ cos cos A 


X = 


Xj 


= 


cos <f> sin A 




Xk 




sin ip 



_ x* 

B_ KT 



where 



sin <f> cos A 
sin <f> sin A 

C08 <f> 



dx 



and 



e = 



XA 



sin A 
cos A 
0 



( 2 ) 



and 



_ <9x 
XA = dX' 



The vectors x, n and e provide the basis for a right-handed orthogonal coordinate system. 



1 



z 




n 




Figure 2. The Course Vector 

Let a course a be designated at a point on the earth’s surface. We wish to determine 
the unit course vector c in the direction a which lies in the tangent plane at the point 
in terms of <f> , A and a (see Fig. 2). An arbitrary vector a, defined at the point P with 
coordinates (r, X t (f>) and lying in the tangent plane at the point, can be rotated clockwise 
(looking from P toward the origin) through an angle a by using the operation Rp(a)a, 
where R p (a) is the rotation operator. R p (a) is a composite rotation in 3-space and can be 
decomposed into fundamental rotations in one of several ways. One way is to proceed as 
follows: First, move P (carrying a with P) along its parallel of latitude to the x-z plane (the 
Greenwich meridian); this is equivalent to a clockwise rotation about the 2 -axis through 
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an angle A and is denoted by R,(A). Next, move P along the Greenwich meridian to 
the equator: this is equivalent to a counterclockwise rotation about the y-axis through an 
angle <f> and is denoted by R y (-<^). The point P now lies on the x-axis with coordinates 
(r, 0, 0) and the vector a at P makes the same angle with the Greenwich meridian as it did 
with respect to the original meridian at (r, A, <f>). Next, leave P on the Greenwich meridian 
at the equator and rotate a through an angle a about the x-axis; this clockwise rotation 
is denoted by R*(a). The vector a has now been rotated through the desired angle a 
with respect to the Greenwich meridian. When P is returned to its original position by 
reversing the steps which got it to the equator on the Greenwich meridian, it will have 
been rotated through the angle a with respect to the original meridian of P ) t.e. R y {<i>) 
followed by R z (-A). The composite rotation operator Rp (ot) is 



Hp(a) = R 3 (-A)R y (0)R*(a)R,(-*)R 3 (A). 



The course vector can then be written as c = Rp(a)n. The fundamental x-, y- and x-axis 
rotation operators are given by 





*1 0 


0 


H.(«) = 


0 cos# sin0 


0 -sin 


9 cos 6 




"cos 9 0 


-sin0" 




0 1 


0 




sin 9 0 


COB 9 




cos 9 


sin0 0" 


R.(«) = 


-sin 6 


cos 9 0 


0 


0 1 



These rotation operators are consistent with a right-handed coordinate system and positive 
signs of 9 for a counterclockwise rotation of the coordinate system or a clockwise rotation 
of a point P about the x-, y- or z-axes, respectively, as viewed looking toward the origin 
from the positive end of the rotation axis [Ref. 3, pg. 43 and Ref. 4, pg. 100]. Some 
simplification in determining Rp(a)n can be obtained by noting that 



R»(”^)Rz(A)n 



0 

0 

1 



= k. 



Thus c can be found from 



c = K,(-A)H,(*)E,(or)k. 



(3) 



If an object’s course vector c is known at some point with position vector p then it is 
easily shown (see Fig. 3) that 

n • c = cos a and 
e • c = sin a. 

EVom these relations, the course a is found to be 

a = qatn(e • c,n • c) (4) 

where qatn is a quadrant determining arctangent function (see Appendix A and Ref. 5). 



n 




Figure 3. An Object’s Course 

The course can also be determined from the great circle vector g that is defined by 

g = x x c, (5) 

where p = rx. From Figure 3 we see that the following relationships hold: 

n • g = cos £ = sin a and 
e g = cos i) = — cos a, 

whence 

a = qatn(n g, -e • g) (6) 

If the object is traveling with speed v and is not maneuvering, it’s course will be a 
great circle. Let vo = vco denote the object’s velocity vector, where Co is its course and 
po is its position vector at time 0. At some subsequent time t, the object’s position vector 
will be p(t) and 

p (f) = op 0 + bv 0 t (7) 
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where a and b are to be determined (see Fig. 4). Dotting Equ. 7 from the right by p 0 , we 
see that 

p(i) • p 0 = ap 0 • p 0 



or, with angles in radians, 



a _ P(<) • Po 
Po • Po 



ip(() ■ Po = ^ 



cos 



(?)] 



or 




Similarly, dotting Equ. 7 from the right by Vo we find 



p(f) • v 0 = &v 0 ■ v 0 f = bv*t 



or 



J =^i p (‘)- T o = 3J 



v 2 t 





so that 
Thus, 




P (0 = Po cos 
In terms of the unit vectors, this equation becomes 



+ v 0 - sin 

v 




rx{t) = rxo cos + vc 0 £ sin ^y^ , 
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or 



x(t) = Xo cos l^pj + Co sin 



( 8 ) 



Applications of the these relations are made in the following sections of this report. 



m. GREAT CIRCLE NAVIGATION. 

The “Direct Solution Algorithm* computes the latitude and longitude of a position P 2 
and the backward azimuth from P 2 to Pj given the latitude and longitude of a position 
Pi, the forward azimuth from Pi to P 2 and the distance from Pi to P 2 . The “Inverse 
Solution Algorithm* computes the distance from position Pi to position P 2 , the forward 
azimuth from Pi to P 2 , and the backward azimuth from P 2 to Pi given the latitude and 
longitude of positions Pi and P 2 . Details of these algorithms using the concept of the 
general spherical triangle are presented in Reference 5. These models are redeveloped here 
using the concepts of Section II. 



A. The Direct Solution Algorithm. Given Ai), forward azimuth (bearing) ai 2 
and distance d , find <^> 2 and A 2 of P 2 and the backward azimuth ar 2 i . Proceed as follows: 

1. Prom <f > i and Ai, compute the components of X] using Equ. 1. 

2. Compute Ci from 

Ci = Ra(~Ai)Ry(0i)R*(oti 2 )k. 



3. Compute x 2 using 

. fd\ 

x 2 = xi cos I - 1 + ci sin I - 1 . 

Note, with d = vt in nautical miles and r = 60(180°/ir), Equ. 8 becomes 



x 2 = Xi cos 




+ Ci sin 




> 



where the arguments of the cosine and sine are now in degrees. 

4. Prom the components of x 2 compute 



fa = sin 1 (xk 2 ) and 
A 2 = qatn(xj 2 ,xj 2 ). 



5. Compute g = Xi xci. 

6. Compute n 2 and e 2 using Equ. 2. 

7. Using Equ. 6 compute a 2i = -qatn(n 2 • g, -e 2 • g). Note that the sign change occurs 
because a 2i is the backward azimuth. 
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B. The Inverse Solution Algorithm. Given Ai) and P 2 (<£ 2 ,A 2 ), find the distance 
d from Pi to P 2 , the forward azimuth ai 2 from Pi to P 2 and the backward azimuth a-n 
from Pi to Pi . Proceed as follows: 

1. From fa and A,- compute x,-, for i = 1, 2. 

2. Compute d = r cos -1 (xi • x 2 ), whence d = 60 ( 180 /tt) cos _1 (xi • x 2 ) is the distance in 
nautical miles. 

3. Compute n, a and e,- for i = 1, 2 (see Equ. 2). 

4. Compute 

_ Xi x Xi 
8 ~ jxi xx 2 |' 

This equation arises because the great circle passes through both Xi and x 2 , hence g 
must be perpendicular to Xi and to x 2 . 

5. Compute ai 2 = qatn(ni • g, -ei • g). 

6. Compute a 2 i = - qatn(n 2 • g, -e 2 • g). 



IV. CLOSEST POINT OF APPROACH (CPA) PROBLEM. 

Consider two objects traveling on different great circle paths. From Equ. 8, their tracks 
will be characterized by the equations 

x,(f) = Xio cos (Jit + Cio sin w,f , for i = 1, 2, (9) 



where w,- = Vi/r. At any time f, their angular separation a(f) in radians is determined by 



coss(f) =X!(t) x 2 (f). 



(10) 



The time of minimum separation (CPA) is that time T when ^(cos^f)] = 0. That 
is, we must find T such that 



m dx 2 (T) dxi{T) 

Xl(T) ' sr + ~dr 



■xi(T) = 0. 



Unfortunately, this equation cannot be solved explicitly. One approach is to use the 
Newton- Raphson iteration method (Ref. 6] to find the root T of /(f), where 



/(f)=Xi(f)- 



dx 2 (f) dxi(f) 
dt dt 



x 2 (f). 



Taking the required derivatives of Equ. 9 and performing the required dot products, we 
find that 



f(t) — A sin oji t sin w 2 f + B cos wi f cos w 2 f + C sin Wi t cos w 2 f + D cos wi f sin a/ 2 f , 
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where 

A = — (wiXjo • C20 4 - U2X20 • C10), 

B = ct^iCio • X20 + W2C20 • 3Cio» 

C = -(wiXio • X20 - W2C20 • C10) and 
D = (J1C10 • C20 — W2X20 • X10. 

The use of the Newton-Raphson method requires that the derivative of f(t) with respect 
to t, namely /'(<), be computed or estimated. We find that 

/'(t) = - ( Cu 2 4 - D<jJi)Bmu)itBm(jj<2t 4 - (Du? 4 - Cwi) coswjt cosa^t 
+ (Aw 2 - 5 wi)smwi<cosw 2 f - (Bw 2 - Aw* ) cos wi* sin 0/2*. 

The Newton-Raphson method requires us to compute 

u * l = u ~7tfi fori = 0 ' 1 >- (11) 

where to is an initial estimate of T t until some value of i is found for which f(U)/f'(U) is 
sufficiently close to zero. 

The CPA option in the computer program (Appendix B) will print out the time to 
CPA (from time t = 0 ), the distance between objects at CPA and the bearing from object 
1 from object 2 at CPA. Also printed is the number of iterations required for convergence 
of Equ. 11 to |/(t*)// f (f*)| less than 0.00001 hours. A negative time to CPA indicates that 
CPA has already occurred. 

V. INTERCEPT PROBLEM. 

The intercept problem is divided into two problems, each of which may require an answer. 
In both problems we are given the initial position, course and speed of a target as well 
as the position of an interceptor or launch platform. In the first problem we are given 
the time (or elapsed time) at which intercept is desired and are required to compute the 
vehicle speed needed for an intercept to take place. In the second problem we are given 
the speed of an intercept vehicle and wish to compute the time required for an intercept 
to occur provided an intercept can be made. Provision for both of these options is made 
in the program presented in Appendix B. 

A. Speed Required to Intercept. Inputs are the initial latitude and longitude of an 
interceptor and a target, and the target course and speed. Also input is the time of 
the desired intercept. Outputs are the speed required of the interceptor, the course of 
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the interceptor, the distance traveled to intercept, and the latitude and longitude of the 
intercept. 

From the inputs, use Equ. 1 to compute xio and x 2 o, the position vectors of the 
interceptor and target, respectively, at time t = 0. Denote the time required to intercept 
by t. Compute the target course vector, C 20 , using Equ. 3 and then compute x 2 (t), the 
position of the target at the time of intercept, using Equ. 8. The remainder of the problem 
is solved using the inverse solution algorithm discussed previously. The speed required for 
intercept is given by Vi = d/t, where d is the distance from the initial interceptor position 
to the target position at the time of intercept. 

B* Time Required to Intercept. Inputs are the initial latitude and longitude of an 
interceptor and a target, and the target course and speed. For a given interceptor speed, 
it may or may not be possible to make an intercept. We develop two sub- algorithms. The 
first algorithm is to compute the minimum interceptor speed required to achieve intercept 
and the time required to make such an intercept. The second algorithm queries the user to 
input an interceptor speed. If the speed is at least that required for intercept, then the time 
required to intercept is computed. (If the interceptor speed is greater than the minimum 
required, there are two possible solutions for the time to intercept — the shortest time to 
intercept is computed by the algorithm). Outputs are the minimum required interceptor 
speed, the time to intercept at minimum speed, the course of the interceptor, the distance 
traveled to intercept, and the latitude and longitude of the intercept. 

The first problem is to determine the minimum speed, v m , required to make an in- 
terception. This can be accomplished by finding the time of intercept t m which makes 
dv/dt = 0. We can relate v to s, the angular separation in radians, between the two points 
Xio and x 2 (t) by the relation t;(f) = rs(t)/t. We find that dv(t)/dt = 0 implies 



r 



‘Ids s 
Jit ~ t? 



= 0 . 



(12) 



Anticipating that it will not be possible to find a closed form solution, we prepare to use the 
Newton- Raphson procedure (Equ. 11). Multiplying Equ. (12) by gives us the function 






for which we wish to find t m such that f[t m ) = 0. Also needed is 
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Using Equ. 10 to determine s(t) we find that 



ds _ 1 dx 3 (f) 

dt * sins * 10 dt 



and 



dh = 1 _ 

dt 2 sin a 



CO 83 
. i 

sm s 




} 



where 



Xio • — ^ — — - u>2 [xio ■ *20 sino^f — Xio ■ ®20 coswat] and 

d 2 xi(t) 21 4 , 

XiQ = -w 2 [Xio • X 20 cosw 2 < 4 - Xi 0 • c 2 o sm Uit\ . 



The Newton- Raphson procedure continues until t m is found, then v m = sr/t m . 

The second problem is to find the time of interception tj when the interceptor’s speed 
Vi is given. Once more the Newton-Raphson procedure is used. Ab before, we can relate 
Vi to a(<), the angular separation in radians, between two points Xio and X2(f) by the 
relation v(t) = r$(t)/t , which tells us that we must require wi = a(f)/f. That is, we wish 
to find ti for which a(fj)/fj equals the constant wi, or equivalently, we wish to find f/ such 
that /(fj) = 0 where 

Also needed is 

The equation for dajdt is the same as that given in the previous paragraph. The remaining 
output is found using the inverse solution algorithm for the points Xio and X2(fj). 
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VI. SAMPLE PROBLEMS 

Master Menu. The master menu for algorithm demonstration program is shown below. 



ALGORITHM DEMO 

1) DIRECT SOLUTION 

2) INVERSE SOLUTION 

3) FIND CPA 

4) SPEED NEEDED TO INTERCEPT 

5) TIME NEEDED TO INTERCEPT 

6) QUIT 



Problem 1 . Suppose you are at San Francisco (latitude 37° 47' north and longitude 
122°25' west), that your initial course is 260° and that you travel a distance of 4000 n. mi. 
What is your final position? Select Option 1 from the master menu. 

DIRECT SOLUTION 

1st LATITUDE DD.MMSS (-S) ? 37.47 
1st LONGITUDE DDD.MMSS (-E) ? 122.25 
INITIAL COURSE DDD.MMSS ? 260 
DISTANCE (n. mi.) ? 4000 

SPHERICAL EARTH DIRECT SOLUTION 
2nd LATITUDE 6° 4 1.9' 

2nd LONGITUDE -172°00.7' 

BACK AZIMUTH 51°35.9' 

PRESS ANY KEY TO CONTINUE 
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Problem 2. Suppose you are at San Francisco (latitude 37° 47' north and longitude 
122° 25' west) and that your destination is Sydney, Australia (latitude 33° 51' south and 
longitude 151° 13' east). How far do you travel, what is your initial course, and what is the 
backward azimuth from Sydney to San Francisco? Select Option 2 from the master menu. 

INVERSE SOLUTION 

1st LATITUDE DD.MMSS (-S) ? 37.47 
1st LONGITUDE DDD.MMSS (-E) 7 122.25 
2nd LATITUDE DD.MMSS (-S) ? -33.51 
2nd LONGITUDE DDD.MMSS (-E) ? -151.13 

SPHERICAL INVERSE SOLUTION 

DISTANCE 6446.3 n.mi. 

FORWARD COURSE 240° 18. 9' 

BACK COURSE 55°45.9' 

PRESS ANY KEY TO CONTINUE 
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Problem 3. Suppose an observer is at 20° north, 60° west traveling on a course of 
010° at a speed of 15 knots. A target is reported to be at 34° north, 50° west on a course 
of 220° at a speed of 300 knots. Assuming that neither observer or target changes course 
or speed, how much time will elapse until CPA and at CPA where will the target be with 
respect to the observer? Select Option 3 from the master menu. 

FIND CPA 

1st LATITUDE DD.MMSS (-S) ? 20 
1st LONGITUDE DDD.MMSS (-E) ? 60 
INITIAL COURSE DDD.MMSS ? 10 
SPEED (knots) ? IS 

2nd LATITUDE DD.MMSS (-S) ? 34 
2nd LONGITUDE DDD.MMSS (-E) ? 50 
INITIAL COURSE DDD.MMSS ? 220 
SPEED (knots) ? 300 

TIME TO CPA - 3h09m48s 
DISTANCE AT CPA = 67.03 n. mi. 

BEARING AT CPA - 304°06.3 / 

NO. ITERATIONS = 3 

PRESS ANY KEY TO CONTINUE 
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As an additional output, a table of observer positions and target positions is given at 
CPA and six equally spaced times before and after CPA. 

FIND CPA 



TIME 


LAT 1 


LONG 1 


DISTANCE 


BEARING (l->2) 


00s 


20°00.0 / 


60° 00. O' 


994.34 


30° 18. 8' 


31m38s 


20° 07. 8' 


59° 58. 5' 


829.45 


29°31.6' 


Ih03ml6s 


20°15.6 / 


59° 57.1' 


664.78 


28° 21. 6' 


Ih34m54s 


20°23.4 / 


59° 56. 6' 


500.56 


26°26 . 3' 


2h06m32s 


20° 31. 2' 


59° 54.1' 


337.43 


22° 39. 8' 


2h38ml0s 


20° 38. 9' 


59° 52. 7' 


178.42 


12° 02. 7' 


3h09m48s 


20°46.7 / 


59°51 .2' 


67.03 


304° 06. 3' 


3h41m27s 


20° 54. S' 


59°49.7' 


178.43 


236° 10.0' 


4hl3m05s 


21° 02. 3' 


59° 48. 2' 


337.43 


225° 33. 2' 


4h44m43s 


21° 10.1' 


59°46.7' 


500.58 


221°47.3' 


5hl6m21s 


21° 17. 9' 


59° 45. 3' 


664.81 


219°62.7' 


6h47mS9s 


21°25.7' 


59° 43. 8' 


829.5 


218°43.7' 


6hl9m37s 


21°33.4' 


69° 42. 3' 


994.41 


217° 67. 6' 



PRESS ANY KEY TO CONTINUE 
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Problem 4. Suppose an observer at 20° north, 60° west wishes to launch an inter- 
ceptor at a target reported to be at 34° north, 50° west on a course of 220° at a speed of 
600 knots. If interception is required to take place 45 minutes (2700 seconds) after launch, 
how fast must the interceptor travel, and where will the intercept take place? (Assume 
that the target does not change course or speed.) Select Option 4 from the master menu. 

SPEED NEEDED TO INTERCEPT (l->2) 

1st LATITUDE DD.MMSS (-S) ? 20 

1st LONGITUDE DDD.MMSS (-E) ? 60 

2nd LATITUDE DD.MMSS (-S) ? 34 

2nd LONGITUDE DDD.MMSS (-E) ? 50 

2nd COURSE DDD.MMSS ? 220 

2nd SPEED (knots) ? 600 

TIME TO INTERCEPT (SECONDS) ? 2700 

SPEED REQUIRED - 730.1 knots 

BEARING TO INTERCEPT = 26°06.9' 

RANGE TO INTERCEPT - 547.5 n. ml. 

LATITUDE OF INTERCEPT = 28° 08. O' 

LONGITUDE OF INTERCEPT - 56°27.6' 

1) CHANGE TIME OF INTERCEPT 

2) NEW PROBLEM 

3) MASTER MENU 
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Problem 5. As in the previous problem, suppose an observer at 20° north, 60° west 
wishes to launch an interceptor at a target reported to be at 34° north, 50° west on a 
course of 220° at a speed of 600 knots. If the maximum speed of the interceptor is 700 
knots, how long will it take before interception can occur, what should be the interceptor’s 
initial great circle heading, and where will the intercept take place? (Assume that the 
target does not change course or speed.) Select Option 5 from the master menu. 



TIME NEEDED TO INTERCEPT (l->2) 



1st LATITUDE DD.MMSS (-S) ? 20 
1st LONGITUDE DDD.MMSS (-E) ? 60 
2nd LATITUDE DD.MMSS (-S) ? 34 
2nd LONGITUDE DDD.MMSS (-E) ? 50 
2nd COURSE DDD.MMSS ? 220 
2nd SPEED (knots) ? 600 



MINIMUM SPEED REQUIRED TO INTERCEPT = 52.6 knots 
TIME REQ’D TO INTERCEPT AT MIN SPEED - Ih39m50s 



INTERCEPTOR SPEED (knots) ? 700 



TIME REQUIRED 
BEARING TO INTERCEPT 
RANGE TO INTERCEPT 
LATITUDE OF INTERCEPT 
LONGITUDE OF INTERCEPT 



=* 46m03s 
= 25° 56.1' 

■ 537.2 n.mi. 
= 27° 59. 6' 

■ 55° 34 . V 



1) CHANGE INTERCEPTOR SPEED 

2) NEW PROBLEM 

3) MASTER MENU 



16 



V. REFERENCES. 



1. Chauvenet, Wm., A Treatise on Plane and Spherical Trigonometry , 7th ed., J. B. 
Lippincott & Co., 1869. 

2. “Computer Aided Stationing Tool (CAST)”, Programmable Description Document, 
11 May 1984, Prepared by MANTECH International Corporation 

3. Mueller, I. I., Spherical and Practical Astronomy as Applied to Geodesy , Frederick 
Unger Publishing Co., 1969. 

4. Goldstein, H., Classical Mechanics , Addison- Wesley Press, Inc., 1950. 

5. Shudde, R. H., “Some Navigation and Almanac Algorithms*, Naval Postgraduate 
School Technical Report NPS55-85-023, September 1985. 

6. Hildebrand, F. B., Introduction to Numerical Analysis , McGraw-Hill Book Company, 
Inc., 1956. 



17 



APPENDIX A: The QATN Function 



This routine is the standard arctangent function corrected for quadrant. The quadrant arc- 
tangent function is occasionally implemented as the ATAN2 function, the ANGLE function 
or the Rectangular-to-Polar function. 

Entering variables are the x- and y-coordinates, X and Y . The exiting variable is 
the angle 0, where — it < 0 < x. Use of the quadrant arctangent function is denoted by 
0 = qatn(F,X). 

1. If X ^ 0, go to step 4. 

2. Set 0 = (*■/ 2) * sgn(F). 

3. Go to step 8. 

4. Set 0 = arctan(y/X). 

5. If X > 0, go to step 8. 

6. Set 0 = 0 + ir * sgn(F). 

7. If Y = 0, set 0 = x. 

8. Return. 

Note: 

If Y > 0 then sgn(F) = +1. 

If Y = 0 then sgn(F) - 0. 

If Y < 0 then sgn(F) = -1. 

Users of Microsoft BASIC can simplfy the qatn function significantly by using the code 
given below. To return an angle of 0 (designated by A in the code) in the range of (-tt, it), 
use: 

PI = 4*ATN(1) : TP = PI + PI: EPS = IE-33 

A =» ATN (Y/ (X— EPS* (X=0) ) ) - PI*(X<0)*(SGN(Y) - (Y=0)) 

To return a value of A in the range of (0, 2ir), use: 

PI =» 4*ATN(1) : TP = PI + PI: EPS = IE-33 

A = ATN (Y/ (X-EPS* (X=0) ) ) - PI*(X<0) + TP*(X >= 0)*(Y<0) 
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APPENDIX B: Program Listing 



10 * *RECT COORD ALGORITHMS* R.H. SBUDDE, 03-03-86. REV. 03-19-86 1000 

13 ' RECTALGR* ,A 

100 DEFDBL A-Z 

110 PI4-ATN(1#) :PI2-PI4+PI4:PI-PI2+PI2:TP-PI*PI:RD-PI/180#:EP8-lD-33 
120 AE=80#*360#/TP 

130 DEF FNM(X)-X-MO*INT(X/MO) : ’ X MOD MO FUNCTION. 

140 DEF FNL(X)-X-TP*INT((X+PI)/TP): ' LONGITUDE ADJUST (-PI. PI) 

160 DEF FNR(X)"INT(X*M0+.6)/M0: * ROUNDING FUNCTION. 

160 DEF FNG(X)=X+PI*SGN(X)*(AB8(X)>PI2) : * LATITUDE ADJUST (-PI/2, PI/2) 

170 DEF FNAC8(X)"ATN(SQR(1#-X*X)/(X-EPS*(X«0#)))-PI*(X<0#) : * ARCCOS 
180 DEF FNASN(X)=ATN(X/(SQR(1#-X*X)-EPS*(AB8(X)=1#))) : * ARC8IN 
190 *QATN (-PI, PI) FUNCTION: 

200 DEF FNATN2(Y,X)=ATN(Y/(X-EPS*(X=0#)))-PI*(X<0#)*(SGN(Y)-(Y=0#)) 

210 *QATN (O.TWOPI) FUNCTION: 

220 DEF FNATNP(Y.X)=ATN(Y/(X-EPS*(X=0#)))-PI*(X<0#)+TP*(X>=0#)*(Y<0#) 

230 * CROSS PRODUCT : X( . )-Xl ( . ) X X2( . ) : 

240 DEF FNCX(X1,Y1,Z1,X2,Y2,Z2)-Y1*Z2-Z1*Y2 
250 DEF FNCY(X1,Y1,Z1,X2,Y2,Z2)-X2*Z1-Z2*X1 
260 DEF FNCZ(X1,Y1,Z1,X2,Y2,Z2)-X1*Y2-Y1*X2 
270 GOTO 2000 
280 * 

1000 * DECIMAL TO HH MM 88 

1010 V|-" *:IF X<0 THEN V$-»-“:X— X 

1020 X=X+1/7200# : Y=INT(X) : Z=LEN(STR| (Y) ) - 1 

1030 KX-0:IF Y<>0 THEN V|-V|+RIGHT|(* *+8TR|(Y) ,Z)+*h“ :KX-1 

1040 X=60i* (X-Y) : Y=INT(X) 

1050 IF Y<>0 OR KX-1 THEN X|-STR|(100#+Y) :V|-V|*RIGHT|(X|,2)+"m" 

1060 X=60#*(X-Y) :Y=INT(X) :X$=8TR$(100#+Y) :V$ a V$+RIGHT$(X$,2)+*s* : RETURN 
1070 * 

1080 ' DECIMAL TO DDD MM.F 

1090 VI-" ":IF X<0 THEN V$-“-*:X— X 

1100 X=X+1/ 1200# : Y»INT(X) : V|»V$+RIGHT$ ( “ *+STR|(Y) ,3)+CHR$(248) 

1110 X»600f*(X-Y):Y-INT(X):X$-STR$(1000+Y) 

1120 V$=V$+MID$ (X$ ,3,2)+* . *+RIGHT$(X$, 1)+** ■ : RETURN 
1130 * 

1140 * DDD.MMSS TO DECIMAL 

1150 IX-0:FOR Zl-l TO LEN(V$) :C$-MID$(V$,Z! ,1) :IF C$-"."THEN IX-ZI 
1160 NEXT: IF IX=0 THEN X-VAL(V$) : RETURN 
1170 X-VAL(LEFT$(V$,n)):SN-l:IF X<0# THEN SN—SN:X—X 
1 180 V$»V$* *0000* : Y=VAL (MIDI (V| . IX* 1 , 2 ) ) : Z»VAL (MIDI ( V| , IX*3 . 2 ) ) 

1190 X-SN*((Z/60#+Y)/60#+X): RETURN 
1200 * 

1300 ‘DIRECT SOLUTION, SPHER EARTH -RECT COORD. ALL ANGLES MUST BE IN RADIUS 
1310 ‘INPUT: LATITUDE PI. LONGITUDE LI. FORWARD AZIMUTH A12 AND 
1320 ‘ DISTANCE DD TO A POINT P2. NOTE: DD IS IN RADIANS. 

1330 ‘OUTPUT: LATITUDE P2, LONGITUDE L2 AND BACKWARD AZIMUTH A21. 

1340 P=P1:L=-L1:G08UB 1760: 'CHNG SIGN OF LI GIVES RIGHT-HANDED COORDS 
1360 FOR I!=l TO 3:P0S1(I! )=POSI(I! ) :N0RTH1(I!)=N0RTH(I!) :EAST1(I! )=EAST(I! ) 
NEXT 

1360 G08UB 1810:CD=COS(DD):SD=8IN(DD) 
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1370 

1330 

1300 

1400 

1410 

1420 

1500 

1510 

1620 

1530 

1540 

1660 

1660 

1670 

1580 

1500 

1600 

1610 

1620 

1630 

1640 

1660 

1700 

1710 

1720 

1730 

1740 

1760 

1760 

1770 

1780 

1700 

1800 

1810 

1820 

1830 

1840 

1860 

1860 

1870 

2000 

2010 

2020 

2030 

2040 

2060 

2060 

2070 

2080 



FOR I!-l TO 3:P0S2(I!)»P0S1(I!)*CD+CVEC1(I!)*SD:NEXT 

L2=FNATN2(P0S2(2) ,P0S2(1)) :P2-FNASN(P0S2(3)) :P-P2:L»L2:G08UB 1750 

X-GVEC(l) *EAST(1) +GVEC(2) *EAST(2)+GVEC(3) *EAST(3) 

Y=- (GVEC(l) *N0RTH(1)+GVEC(2) *N0RTH(2)+GVEC(3) *N0RTH(3) ) 

A21-FNATNP(Y,X):L2— L2: RETURN: 'CONVERT BACK TO LEFT-HANDED OUTPUT 
» 

'INVERSE SOLUTION, SPHER EARTH -RECT COORD. ALL ANGLES MUST BE IN RADIANS 
'INPUT: LATITUDES PI * P2. AND LONGITUDES LI ft L2. 

'OUTPUT: DISTANCE DD TO A POINT P2. (NOTE: 0 <= DD <= PI RADIANS). 

' FORWARD AZIMUTH A12, AND BACKWARD AZIMUTH A21. 

P=P1 : L=-L1 : GOSUB 1760:' CHNG SIGN OF LI GIVES RIGHT-HANDED COORDS 

FOR I!-l TO 3:P0S1(I!)-P0SI(I!):N0RTH1(I!)«N0RTH(I!):EAST1(I!)-EAST(I!): 

NEXT 

P=P2 : L=-L2 : G08UB 1750:' CHNG SIGN OF L2 GIVES RIGHT-HANDED COORDS 

FOR I!-l TO 3 : P0S2 ( I ! )-POSI (II): N0RTH2 (I 1 )-NORTH(I ! ) : EAST2 (I ! )«EAST(I ! ) : 

NEXT 

DD-FNACS(P0S1(1)*P0S2(1)+P081(2)*P082(2)+P081(3)*P0S2(3)) 

X=FNCX(P0S1(1) ,P0S1(2) ,P0S1(3) ,P0S2(1) ,P0S2(2) ,P082(3)) 

Y-FNCY(P0S1 ( 1) . P0S1 (2) , P0S1 (3) . P0S2 ( 1) ,P0S2(2) , P0S2(3) ) 

Z-FNCZ (P0S1 ( 1 ) , P081 (2) , P0S1 (3) , P0S2 ( 1) , P0S2(2) , P082 (3) ) 

A1 2-FNATNP(X* N0RTH1 ( 1)+Y*N0RTH1 (2)+Z*N0RTHl (3) , 

- (X+EAST1 (1 ) +Y*EAST1 (2) +Z*EAST1 (3) ) ) 
A21-FNATNP(-(X*N0RTH2(1)+Y*N0RTH2(2)+Z*N0RTH2(3)) , 

X+EAST2 ( 1 ) +Y*EAST2 (2 ) +Z* EAST2 (3) ) 

RETURN 

• 

CA-COS(A) :SA"SIN(A) :T«Y*CA+Z*SA: Z«Z*CA-Y*SA:Y"T: RETURN: ’ X-AXIS ROT 

CA=COS (A) :SA=SIN(A):T=Z*CA+X*SA:X=X*CA-Z*SA:Z=T: RETURN : ' Y-AXIS ROT 

CA-C0S(A):SA-SIN(A):T-X*CA+Y*SA:Y-Y*CA-X*8A:X«T:RETURN: ' Z-AXIS ROT 
$ 

'UNIT VECTORS: POSITION, NORTH ft EAST. 

SL=SIN(L) :CL=COS(L) : SP=SIN(P) :CP=COS(P) 

POSI( 1) =CP*CL : POSI (2) =CP*SL : POSI (3) =SP 

N0RTH(1)=-SP*CL:N0RTH(2)=-SP*SL:N0RTH(3)=CP 

EAST( 1)—SL : EAST(2) -CL : EAST (3) -0 : RETURN 
« 

'VECTORS :CVEC=COURSE ft GVEC=GREAT CIRCLE NORMAL 

X=0 : Y=0 : Z=1 : A=A12 : GOSUB 1700 :A=P: GOSUB 1710 :A=-L: GOSUB 1720 

CVEC1 ( 1 ) -X : CVEC1 (2) -Y : CVEC1 (3) -Z 

GVEC(1)»FNCX(P0S1(1) ,P0S1(2),P0S1(3) .CVECl(l) ,CVEC1(2) ,CVEC1(3)) 
GVEC(2)-FNCY(P0S1(1) ,P0S1(2) ,P0S1(3) .CVECK1) ,CVEC1(2) ,CVEC1(3)) 
GVEC(3)=FNCZ(P0S1 (1) , POSI (2) , POSI (3) , CVEC1 (1) , CVECK2) , CVEC1(3)) 

RETURN 

» 

CLS: PRINT SPC(20) ; -ALGORITHM DEMO 

PRINT: PRINT: PRINT 

PRINT SPC(16);-1) DIRECT SOLUTION 

PRINT SPC(16) ; -2) INVERSE 80LUTI0N 

PRINT SPC(16);-3) FIND CPA 

PRINT SPC(16);-4) SPEED NEEDED TO INTERCEPT 

PRINT SPC(16);-6) TIME NEEDED TO INTERCEPT 

PRINT : PRINT SPC(16);-6) QUIT 

GOSUB 9010:C-VAL(C$):0N C GOSUB 3000,4000,6000,6000,7000,8000 
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2090 

2100 

3000 

3010 

3020 

3030 

3040 

3060 

3060 

3070 

3080 

3090 

3100 

3110 

3120 

3130 

3140 

3160 

4000 

4010 

4020 

4030 

4040 

4060 

4060 

4070 

4080 

4090 

4100 

4110 

4120 

4130 

4140 

4160 

4160 

4170 

4180 

6000 

6010 

6020 

6030 

5040 

5050 

6060 

6070 

6080 

5090 

5100 

6110 

5120 

5130 

6140 

5150 



GOTO 2000 



CLS: PRINT 8PC(15) ; "DIRECT SOLUTION" : PRINT : PRINT 
PRINT SPC(IO); : PRINT "lat LATITUDE DD.MM88 (-8) ■; 

INPUT V$:GOSUB 1150:P1-X*RD 

PRINT SPC(IO); : PRINT "lat LONGITUDE DDD.MMSS (-E) "; 

INPUT V$ : GO SUB 1150.L1-X+RD 

PRINT SPC(IO); : PRINT "INITIAL COURSE DDD.MMSS ■; 

INPUT V$ : GO SUB 1150:A12»X*RD 

PRINT SPC(IO);: INPUT "DISTANCE (n. mi.) ? ",D 

Dl-D+RD/60# 

M0-100 

DD"D1 : GO SUB 1340: PRINT: PRINT SPC (8) ; "SPHERICAL EARTH DIRECT SOLUTION 
PRINT SPC(12) ; "2nd LATITUDE "; :X»P2/RD: GOSUB 1090:PRINT V$ 

PRINT SPC(12) ;"2nd LONGITUDE ■ ; :X-L2/RD:G0SUB 1090:PRINT V| 

PRINT SPC(12);"BACK AZIMUTH "; :X-A21/RD:G0SUB 1090:PRINT V| 

GOSUB 9000: GOTO 2000 

a 

CLS: PRINT SPC (1 5) ; "INVERSE SOLUTION" : PRINT: PRINT 
PRINT SPC(IO); : PRINT "lat LATITUDE DD.MMSS (-S) "; 

INPUT V$: GOSUB 1150:P1-X 

PRINT SPC(IO); : PRINT “lat LONGITUDE DDD.MMSS (-E) "; 

INPUT V$: GOSUB 1150:L1-X 

PRINT SPC(IO); : PRINT "2nd LATITUDE DD.MMSS (-S) "; 

INPUT VI: GOSUB 1150:P2-X 

PRINT SPC(IO); : PRINT "2nd LONGITUDE DDD.MMSS (-E) ■; 

INPUT V$ :GOSUB 1150.L2-X 

Pl-Pl *RD : P2=P2*RD : Ll-Ll *RD : L2»L2*RD 

DD-D1: GOSUB 1640 

PRINT: PRINT SPC (8) ; "SPHERICAL INVERSE SOLUTION 
• 

PRINT SPC(12); "DISTANCE "; 

M0*100: PRINT FNR(60#*DD/RD) ; • n.mi. 

PRINT SPC(12); "FORWARD COURSE ■ ; :X-A12/RD:G0SUB 1090:PRINT V$ 

PRINT SPC(12);"BACK COURSE ■ ; :X»A21/RD: GOSUB 1090:PRINT V$ 

GOSUB 9000: GOTO 2000 
» 

CLS: PRINT SPC (20) ; "FIND CPA" : PRINT : PRINT 

PRINT SPC(IO); : PRINT "lat LATITUDE DD.MMSS (-S) •; 

INPUT V$: GOSUB 1160:P1-X*RD 

PRINT SPC(IO);: PRINT "lat LONGITUDE DDD.MMSS (-E) "; 

INPUT V$: GOSUB 1150:L1»X*RD 

PRINT SPC(IO); .PRINT "INITIAL COURSE DDD.MM8S ■; 

INPUT V$: GOSUB 1150:A1»X*RD 

PRINT SPC(IO);: INPUT "SPEED (knota) ? ",S1 

PRINT SPC(IO); -.PRINT "2nd LATITUDE DD.MMSS (-S) ■; 

INPUT V$: GOSUB 1150:P2=X*RD 

PRINT 8PC(10); : PRINT "2nd LONGITUDE DDD.MMSS (-E) "; 

INPUT V$: GOSUB 1150:L2»X*RD 

PRINT SPC(IO); : PRINT "INITIAL COURSE DDD.MMSS "; 

INPUT V$: GOSUB 1160:A2-X*RD 

PRINT SPC(IO); : INPUT "SPEED (knota) ? ",S2 

B1-S1/AE:B2-S2/AE 
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5160 

5170 

5180 

5190 

5200 

5210 

5220 

5230 

5240 

5250 

5260 

5270 

5280 

5290 

5300 

5310 

5320 

5330 

5340 

5350 

5360 

5370 

5380 

5390 

5400 

5410 

5420 

5430 

5440 

5460 

5460 

5470 

5480 

5490 

5600 

5610 

5520 

5530 

5540 

5550 

5660 

5570 

5580 

5590 

5600 

5610 

5620 

5630 

5640 

5650 

5660 

5670 



P=P1 : L=-L1 : GOSUB 1750:* CHNG 8IGN OF LI GIVES RIGHT-HANDED COORDS 
FOR I!-l TO 3:X1(I!)-P0SI(I!):NEXT 

X=0 : Y=0 : Z=1 : A=A1 : GOSUB 1700 :A=P: GOSUB 1710 :A=-L: GOSUB 1720 
C1(1)-X:C1(2)-Y:C1(3)-Z 

P=P2:L=-L2: GOSUB 1750:' CHNG SIGN OF L2 GIVES RIGHT-HANDED COORDS 
FOR I!-l TO 3:X2(I!)-P0SI(I!):NEXT 

X»0:Y-0:Z-1:A-A2: GOSUB 1700.A-P: GOSUB 1710 : A— L: GOSUB 1720 

C2(1)-X:C2(2)-Y:C2(3)-Z 

• 

X1C2-X1 ( 1 ) *C2 ( 1 ) +X1 (2) *C2 (2) +X1 (3) *C2 (3) 

C1X2=C1 (1)*X2(1)+C1(2)*X2(2)+C1(3)*X2(3) 
X1X2-X1(1)*X2(1)+X1(2)*X2(2)+X1(3)*X2(3) 
C1C2=C1(1)*C2(1)+C1(2)*C2(2)+C1(3)*C2(3) 

BA=-B1*X1C2 - B2*C1X2 
BB- B1+C1X2 + B2*X1C2 
BO-B1+X1X2 + B2+C1C2 
BD- B1*C1C2 - B2*X1X2 

I 

T=1 : IT ! =0 : ‘ ITERATE WITH NEWT0N-RAPH80N 

B1T-B1*T:B2T-B2*T:S1-SIN(B1T) :C1-C0S(B1T) :S2-SIN(B2T) :C2-C08(B2T) 

S1S2»S1*S2:C1C2*C1*C2:S1C2=S1*C2:C1S2«C1*S2 

F=BA*S1S2+BB*C1C2+BC*81C2+BD*C1S2 

FP— (BC*B2+BD*B1) *8182+(BD*B2+BC*B1) *C1C2+(BA*B2-BB*B1) *S1C2- 
(BB*B2-BA*B1) *C1S2 

IT ! "IT ! ♦ 1 : CORR"F/FP : T»T-CORR : IF AB8(C0RR)<. 00001 THEN 5440 

IF IT1>50 THEN PRINT “NO CONVERGENCE" : STOP 

GOTO 5360 
» 

B1T-B1*T:B2T-B2*T:CB1T-C0S(B1T) :SB1T-SIN(B1T) :CB2T-C0S(B2T) :SB2T-8IN(B2T) 
FOR I!-l TO 3 : P0S1 (I ! )"X1 (I ! ) +CB1T+C1 (I ! ) *SB1T 
P082 ( I ! ) -X2 ( I ! ) *CB2T+C2 ( I ! ) +SB2T : NEXT I ! 

P1=FNASN(P0S1 (3) ) : L1=-FNATN2(P081 (2) , P081 (1) ) 

P2-FNASN(P0S2(3)) :L2— FNATN2(P0S2(2) ,P0S2(1)) 

GOSUB 1540 

X b T: GOSUB 1010: PRINT: PRINT SPC(IO) ; “TINE TO CPA ■ “;V$ 

M0=100#: PRINT SPC( 10) ; “DISTANCE AT CPA » “ ;FNR(60#*DD/RD) ; “ n.al. 

PRINT 8PC(10); “BEARING AT CPA - “; :X-A12/RD: GOSUB 1090:PRINT V$ 

PRINT SPC(IO) ; “NO. ITERATIONS - “ ; IT ! 

TCPA=T: GOSUB 9000 

l 

CLS: PRINT 8PC(22) ;"FIND CPA“:PRINT 

PRINT “ TIME LAT 1 LONG 1 DISTANCE BEARING(l->2) “ 

DT"TCPA/6# : T"0 : FOR T!-l TO 13:B1T-B1*T:B2T-B2*T 
FOR I!-l TO 3:P081(I!)-X1(I!)*C08(B1T)+C1(I!)*SIN(B1T) 

P0S2 ( I ! )"X2(I ! ) *C03(B2T) +C2(I ! ) *8IN(B2T) : NEXT I! 
P1=FNA8N(P081(3)):L1=-FNATN2(P0S1(2),P0S1(1)) 

P2-FNASN(P082(3)) : L2— FNATN2(P0S2(2) , P0S2(1)) 

GOSUB 1540 

X-T: GOSUB 1010 -.PRINT V$; 

LOCATE CSRLIN, 12 :X-P1/RD: GOSUB 1090:PRINT V$; 

LOCATE CSRLIN, 23 :X-L1/RD: GOSUB 1090.PRINT V$; 

LOCATE CSRLIN ,37: PRINT FNR(60#*DD/RD) ; 
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6680 

6690 

6700 

6710 

6000 

6010 

6020 

6030 

6040 

6060 

6060 

6070 

6080 

6090 

6100 

6110 

6120 

6130 

6140 

6160 

6160 

6170 

6180 

6190 

6200 

6210 

6220 

6230 

6240 

6260 

6260 

6270 

6280 

6600 

6610 

6620 

6630 

6640 

6660 

6660 

6670 

6680 

6690 

6600 

6610 

6620 

6630 

6800 

6810 

6820 

6830 

6840 

6860 



LOCATE CSRLIN , 62 : X"A12/RD : GOSUB 1090: PRINT V$ 

T*T+DT : NEXT T! 

GOSUB 9000: GOTO 2000 
• 

CLS: PRINT SPC ( 10) ; -SPEED NEEDED TO INTERCEPT (l->2)“ :PRINT:PRINT 
CLNX-0 : GOSUB 6610: ’ GET INPUT 

CLNX- 10 -.LOCATE CLNX,1:F0R 1 1-1 TO 11:PRINT SPC(79) :NEXT: ’ CLEAR SCREEN 

LOCATE CLNX, 11 :PRINT -TIME TO INTERCEPT (SECONDS) -;:INPUT TNI 

TMI-TMI/3600# : TM-TMI : PI-PIS : L1-L1S : P2-P2S : L2-L2S : A2-A2S 
• 

• COMPUTE SPEED 

P=P1S:L=-L1S: GOSUB 1760:* CHNG SIGN OF LI GIVES RIGHT-HANDED C00RD8 
FOR 1 1 — 1 TO 3:X1(I!)-P0SI(II):NEXT 

P»P2S:L— L2S: GOSUB 1760: • CHNG SIGN OF L2 GIVES RIGHT-HANDED COORDS 
FOR Ii-1 TO 3 : X2(I I )— POSICI I ) :NEXT 

X=0 : Y=0 : Z=1 : A=A2 : GOSUB 1700 :A=P: GOSUB 1710 :A»-L: GOSUB 1720 

C2(1)»X:C2(2)»Y:C2(3)*Z 

B2T=B2*TN:CB2T=C0S(B2T) :SB2T=SIN(B2T) :CS=0# 

FOR I 1-1 TO 3:P0S2(I!)-X2(I!)*CB2T*C2(II)*SB2T 
C8*CS+X1(I! )*P0S2(I1 ) :NEXT I! 

S-FNACS (CS) : SPD-8+AE/ (TM-EPS* (TM-O) ) 

f 

• GET INVERSE 80LN 

P2-FNASN(P0S2(3)) :L2— FNATN2(P082(2) ,P0S2(1)) : GOSUB 1640 

I 

M0-10: LOCATE CLNX+2, 11 : PRINT -SPEED REQUIRED ■ -;FNR(SPD);“ knota" 
GOSUB 6810: * PRINT OUT BEARING, RANGE, LAT k LONG 
LOCATE CLNK+8, 16: PRINT -1) CHANGE TIME OF INTERCEPT" 

LOCATE CLNX+9. 16 -.PRINT -2) NEW PROBLEM- 
LOCATE CLNX+10,16:PRINT -3) MASTER MENU- 
COSUB 9010:C 3 VAL(C$) : ON C GOTO 6020,6000,2000 

GOTO 6260 

• 

• INPUT ROUTINE 

LOCATE CLN%+3, 11: PRINT "lat LATITUDE DD.MMSS (-S) 

INPUT V$: GOSUB 1160:P1S=X*RD 

LOCATE CLNX+4, 11: PRINT -lat LONGITUDE DDD.MMSS (-E) -; 

INPUT V|: GOSUB 1150:L1S-X*RD 

LOCATE CLNX+6, 11: PRINT -2nd LATITUDE DD.MMSS (-S) 

INPUT V$: GOSUB 1160:P2S=X*RD 

LOCATE CLNX+6, 11: PRINT “2nd LONGITUDE DDD.MMSS (-E) -; 

INPUT V$: GOSUB 1160:L2S»X*RD 

LOCATE CLNX+7, 11: PRINT “2nd COURSE DDD.MMSS •; 

INPUT V$: GOSUB 1160:A2S-X*RD 

LOCATE CLNX+8, 11: INPUT -2nd SPEED (knota) ? “,S2 

B2-S2/AE: RETURN 

• 

* OUTPUT BEARING, RANGE, LAT k LONG 

LOCATE CLNX+3, 11: PRINT -BEARING TO INTERCEPT = -; 

X=A12/RD: GOSUB 1090: PRINT V$ 

LOCATE CLNX+4. 11: PRINT -RANGE TO INTERCEPT = -;FNR(60#*DD/RD) ; - n.mi." 
LOCATE CLNX+5, 11: PRINT -LATITUDE OF INTERCEPT = -; 

X=P2/RD: GOSUB 1090: PRINT V$ 
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6860 

6870 

6880 

6800 

7000 

7010 

7020 

7030 

7040 

7060 

7060 

7070 

7080 

7090 

7100 

7110 

7120 

7130 

7140 

7160 

7160 

7170 

7180 

7190 

7200 

7210 

7220 

7230 

7240 

7260 

7260 

7270 

7280 

7290 

7300 

7310 

7320 

7330 

7340 

7360 

7360 

7370 

7380 

7390 

7400 

7410 

7420 

7430 

7440 

7460 

7460 

7470 

7480 



LOCATE CLNX+6, 11: PRINT "LONGITUDE OF INTERCEPT - ■; 

X-L2/RD: GOSUB 1090: PRINT V$ 

RETURN 

$ 

CL8: PRINT 8PC(10) ; "TIME NEEDED TO INTERCEPT (l->2)": PRINT: PRINT 

CLN%=0 : GOSUB 6610:’ GET INPUT 
* 

* FIND Vain AND Tain vel. 

* Compute arc distance 8 

P=P1S:L=-L1S: GOSUB 1760:’ CHNG SIGN OF LI GIVES RIGHT-HANDED COORDS 
FOR I!=l TO 3:X1(I!)=P08I(I!):NEXT 

P-P28 : L=-L2S : GOSUB 1760:' CHNG SIGN OF L2 GIVES RIGHT-HANDED COORDS 
FOR I!=l TO 3:X2(I!)»P0SI(I!):NEXT 

X=0 : Y=0 : Z=1 : A=A28 : GOSUB 1700 :A=P: GOSUB 1710 : A— L: GOSUB 1720 
C2(1)-X:C2(2)-Y:C2(3)-Z 

f 

' INITIALIZE 

X1X2»X1(1)*X2(1)+X1(2)*X2(2)+X1(3)*X2(3) 

X1C2=X1(1)*C2(1)+X1(2)*C2(2)+X1(3)*C2(3) 

9 

’ BEGIN ITERATION 

S-FNAC8 (X1X2) : T-S*AE/S2 : IT ! -1 

B2T-B2+T : CB2T«C08(B2T) : 8B2T-SIN(B2T) 

CS-X1X2*CB2T+X1C2*SB2T : S-FNAC8 (CS) : SS-SIN (S) 

DSDT=(X1X2*8B2T-X1C2*CB2T) *B2/SS 
F=T*DSDT-S 

FP-T*(B2*B2*(X1X2*CB2T+X1C2*SB2T)-CS*DSDT*DSDT)/SS 
IT!=IT!+1 :CORR=F/FP:T=T-CORR:IF ABS(CORR)<. 00001 THEN 7270 
IF IT! >60 THEN PRINT "NO CONVERGENCE" : STOP 
GOTO 7180 

9 

TMS»T:VMIN«S*AE/T 
LOCATE CLNl+10, 11 : M0-10 

PRINT "MINIMUM SPEED REQUIRED TO INTERCEPT - -;FNR(VMIN) ;■ knots" 

LOCATE CLNX+11, 11: X-TMS: GOSUB 1010 

PRINT "TIME REQ'D TO INTERCEPT AT MIN SPEED - ";V$ 

9 

CLNX-13: LOCATE CLNX,1:F0R I!-l TO 11: PRINT SPC(79) :NEXT: ’ CLEAR SCREEN 
LOCATE CLNX,11:PRINT "INTERCEPTOR SPEED (knots) ";:INPUT 8PD 
IF SPD>=VMIN THEN 7390 

LOCATE CLNX+2, 11: PRINT "SPEED TOO LOW, CANNOT INTERCEPT" 

GOSUB 9000: GOTO 7330 

9 

Bl-SPD/AE:T-TMS/5: T-.l :IT!-1 

B2T*B2*T : CB2T-C0S (B2T) : SB2T»8IN (B2T) 

CS-X1X2*CB2T+X1C2*SB2T:S»FNACS(CS) : SS-SIN (S) 

D8DT=(X1X2*8B2T-X1C2*CB2T) *B2/8S 
F-S/T-Bl :FP-(DSDT-8/T)/T 

IT ! -IT ! + 1 : CORR=F/FP : T-T-CORR : IF ABS(CORR)<. 00001 THEN 7500 
IF IT! >60 THEN PRINT "NO CONVERGENCE" : STOP 
IF ABS(CORR)< 1000000000# THEN 7400 

LOCATE CLNX+2, 11: PRINT "SPEED TOO HIGH, NO CONVERGENCE" 

LOCATE CLNX+4 : GOSUB 9000: GOTO 7330 
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7400 

7500 

7510 

7520 

7530 

7540 

7550 

7560 

7570 

7580 

7590 

7600 

7610 

7620 

7630 

8000 

8010 

9000 

9010 

9020 

9030 



X=T:G08UB 1010 

LOCATE CLNV2, 11: PRINT "TINE REQUIRED - ";V$ 

a 

1 GET INVERSE SOLN 

FOR I!-l TO 3 : P082 ( I ! )*X2 (II) +CB2T+C2 ( I ! ) * SB2T : NEXT I! 

P1=P1S : L1=L1S : P2=FNASN(P0S2 (3) ) : L2=-FNATN2 (P0S2 (2) , P0S2 ( 1) ) : GOSUB 1640 

a 

GOSUB 6810: ' PRINT OUT BEARING, RANGE, LAT & LONG 
LOCATE CLNX+8, 15 .PRINT "1) CHANGE INTERCEPTOR SPEED" 

LOCATE CLNX+9, 15 .PRINT "2) NEW PROBLEM" 

LOCATE CLN%+10, 15: PRINT "3) MASTER MENU" 

GOSUB 9010 : C»VAL(C$) : ON C GOTO 7330,7000,2000 
GOTO 7610 

a 

CLS:END 

a 

PRINT: PRINT SPC( 10) ; "PRESS ANY KEY TO CONTINUE 
FOR I!=l TO 9 : C$= INKEY) : NEXT 
C$=INKEY$ : IF C$="" THEN 9020 
RETURN 
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